Validation

import os
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio
from rasterio import features

import matplotlib.pyplot as plt
import sklearn
from sklearn.metrics import confusion_matrix, f1_score, accuracy_score, classification_report

from pathlib import Path
from IPython.display import display

import plotly.express as px
import plotly.figure_factory as ff
import plotly.graph_objects as go
import plotly.offline
plotly.offline.init_notebook_mode()

print('All libraries successfully imported!')
print(f'Scikit-learn: {sklearn.__version__}')

All libraries successfully imported!
Scikit-learn: 0.24.2

Set directory

computer_path = '/export/miro/ndeffense/LBRAT2104/'
grp_letter    = 'X'

lut_path  = f'{computer_path}data/LUT/'

# Directory for all work files
work_path = f'{computer_path}GROUP_{grp_letter}/WORK/'

in_situ_path = f'{work_path}IN_SITU/'
classif_path = f'{work_path}CLASSIF/'

am_path = f'{work_path}ACCURACY_METRICS/'


Path(am_path).mkdir(parents=True, exist_ok=True)

print(f'Accuracy Metrics path is set to : {am_path}')

Accuracy Metrics path is set to : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/ACCURACY_METRICS/

Set parameters

site = 'NAMUR'
year = '2020'

feat_nb = 2
no_data = 0

ws = 3          # Window size (filtering post classification)

reclassif_flag = False
filter_flag    = False

# Field used for classification
field_classif_code = 'grp_1_nb'
field_classif_name = 'grp_1'

# Field used for reclassification
field_reclassif_code = 'grp_A_nb'
field_reclassif_name = 'grp_A'

Set filenames

s4s_lut_csv = f'{lut_path}crop_dictionary_new.csv'

in_situ_val_shp = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_VAL.shp'
in_situ_val_tif = f'{in_situ_path}{site}_{year}_IN_SITU_ROI_VAL.tif'

if not reclassif_flag and not filter_flag:
    print(f'Reclassification : {reclassif_flag}')
    print(f'Filter           : {filter_flag}')
    classif_tif   = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}.tif'

elif reclassif_flag and not filter_flag:
    print(f'Reclassification : {reclassif_flag}')
    print(f'Filter           : {filter_flag}')
    classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}_reclassify_{field_reclassif_name}.tif'

elif reclassif_flag and filter_flag:
    print(f'Reclassification : {reclassif_flag}')
    print(f'Filter           : {filter_flag}')
    classif_tif = f'{classif_path}{site}_{year}_classif_RF_feat_{feat_nb}_{field_classif_name}_reclassify_{field_reclassif_name}_filter_ws_{ws}.tif'

else:
    print('No classfication file is available !')
    classif_tif = ""

# Confusion matrix and Accuracy metrics files
cm_csv  = f'{am_path}{os.path.basename(classif_tif)}_CM.csv'
cm_html = f'{am_path}{os.path.basename(classif_tif)}_CM.html'
am_html = f'{am_path}{os.path.basename(classif_tif)}_AM.html'


# USED ONLY FOR VISU ON WEBSITE !
#cm_html = f'/export/miro/ndeffense/LBRAT2104/GIT/eo-toolbox/figures/{site}_{year}_CM.html'
#am_html = f'/export/miro/ndeffense/LBRAT2104/GIT/eo-toolbox/figures/{site}_{year}_AM.html'

print(f'Classification file used : \n {classif_tif}')
print(f'Validation polygons used : \n {in_situ_val_shp}')

Reclassification : False
Filter           : False
Classification file used : 
 /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/CLASSIF/NAMUR_2020_classif_RF_feat_2_grp_1.tif
Validation polygons used : 
 /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.shp

1. Prepare data

1.1 Rasterize in situ data validation shapefile

# Open the shapefile with GeoPandas
in_situ_gdf = gpd.read_file(in_situ_val_shp)

# Open the raster file you want to use as a template for rasterize
print(f'Raster template file : {classif_tif}')

src = rasterio.open(classif_tif, "r")

# Update metadata

out_meta = src.meta
out_meta.update(nodata=no_data)

crs_shp = str(in_situ_gdf.crs).split(":",1)[1]
crs_tif = str(src.crs).split(":",1)[1]

print(f'The CRS of in situ data is    : {crs_shp}')
print(f'The CRS of raster template is : {crs_tif}')

if crs_shp == crs_tif:
    print("CRS are the same")

    print(f'Rasterize starts : {in_situ_val_shp}')

    # Burn the features into the raster and write it out

    dst = rasterio.open(in_situ_val_tif, 'w+', **out_meta)
    dst_arr = dst.read(1)

    # This is where we create a generator of geom, value pairs to use in rasterizing

    geom_col = in_situ_gdf.geometry
    code_col = in_situ_gdf[field_classif_code].astype(int)

    shapes = ((geom,value) for geom, value in zip(geom_col, code_col))

    in_situ_arr = features.rasterize(shapes=shapes,
                                     fill=no_data,
                                     out=dst_arr,
                                     transform=dst.transform)
    
    dst.write_band(1, in_situ_arr)

    print(f'Rasterize is done : {in_situ_val_tif}')

    # Close rasterio objects
    src.close()
    dst.close()

else:
    print('CRS are different --> repoject in-situ data shapefile with "to_crs"')

Raster template file : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/CLASSIF/NAMUR_2020_classif_RF_feat_2_grp_1.tif
The CRS of in situ data is    : 32631
The CRS of raster template is : 32631
CRS are the same
Rasterize starts : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.shp
Rasterize is done : /export/miro/ndeffense/LBRAT2104/GROUP_X/WORK/IN_SITU/NAMUR_2020_IN_SITU_ROI_VAL.tif

1.2 Build y_pred and y_true

# Open in-situ used for validation
src = rasterio.open(in_situ_val_tif, "r")
val_arr = src.read(1)
src.close()

# Open classification map
src = rasterio.open(classif_tif, "r")
classif_arr = src.read(1)
src.close()

# Get the postion of validation pixels
idx = np.where(val_arr == no_data, 0, 1).astype(bool)

# Ground truth (correct) target values
y_true = val_arr[idx]

print(f'Reference data (truth) : {y_true}')

# Estimated targets as returned by a classifier.
y_pred = classif_arr[idx]

print(f'Classification data    : {y_pred}')

Reference data (truth) : [   3    3    3 ... 1111 1111 1111]
Classification data    : [   3    3    3 ... 1111 1111 1111]

Sometimes, some classes do not appear in the classification map, they are not predicted by the Random Forest.

This means that some classes in y_true don't appear in y_pred.

classes_all  = sorted(np.unique(y_true))
classes_pred = sorted(np.unique(y_pred))

classes_missing = set(y_true) - set(y_pred)

print(f'{len(classes_missing)} classes are missing in the classification (y_pred) : {classes_missing} \n')

print(f'All training classes :\n {classes_all}')
print(f'All predicted classes (at least once) :\n {classes_pred}')

0 classes are missing in the classification (y_pred) : set() 

All training classes :
 [3, 21, 22, 69, 81, 84, 121, 1111, 1121, 1152, 1171, 1192, 1435, 1511, 1771, 1811, 1911, 1923, 9212]
All predicted classes (at least once) :
 [3, 21, 22, 69, 81, 84, 121, 1111, 1121, 1152, 1171, 1192, 1435, 1511, 1771, 1811, 1911, 1923, 9212]

1.3 Find class names

lut_df = pd.read_csv(s4s_lut_csv, sep=';')

classes_name = lut_df[lut_df[field_classif_code].isin(classes_all)].sort_values(field_classif_code)[field_classif_name].drop_duplicates().to_list()

for code,name in zip(classes_all, classes_name):
    print(f'{code} - {name}')

3 - Grassland and meadows
21 - Fruits trees
22 - Vineyards
69 - Forest
81 - Urban
84 - Greenhouses
121 - Leafy or stem vegetables
1111 - Winter wheat
1121 - Maize
1152 - Barley six-row
1171 - Oats
1192 - Other cereals
1435 - Rapeseed
1511 - Potatoes
1771 - Peas
1811 - Sugar beet
1911 - Alfalfa
1923 - Flax  hemp and other similar crops
9212 - Permanent rivers and streams

2. Confusion Matrix

2.1 Compute Confusion Matrix

cm = confusion_matrix(y_true, y_pred)

cm_df = pd.DataFrame(cm)

cm_values = cm_df.values

cm_df.columns = classes_all
cm_df.index = classes_all

# Export CM in a CSV file
cm_df.to_csv(cm_csv, index=True, sep=',')

display(cm_df)

	3	21	22	69	81	84	121	1111	1121	1152	1171	1192	1435	1511	1771	1811	1911	1923	9212
3	17462	41	17	300	24	0	0	67	0	0	3	2	0	0	0	0	2	0	0
21	686	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0
22	14	40	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
69	48	0	0	732	33	0	0	0	0	0	0	0	0	0	0	0	0	0	0
81	23	3	0	13	404	0	0	1	0	0	0	0	0	0	0	0	0	0	13
84	23	0	0	0	24	2	0	20	0	0	0	0	0	0	0	0	0	0	0
121	5	0	0	0	0	0	64	0	696	0	0	0	0	0	0	20	0	0	0
1111	86	0	7	3	39	1	0	19950	1	1063	47	955	2	1	7	3	0	65	0
1121	63	1	0	0	24	0	1	0	6338	0	2	5	7	11	34	49	0	0	0
1152	256	0	0	0	4	0	0	1182	0	2983	124	424	23	0	0	0	0	0	0
1171	0	0	0	0	0	0	0	12	0	0	90	0	0	0	0	0	0	140	0
1192	17	0	0	1	9	0	0	3960	0	0	112	512	0	0	1	0	0	28	0
1435	76	0	0	0	1	0	0	24	0	86	0	0	463	0	0	0	0	0	0
1511	0	0	3	1	5	0	3	1	130	0	0	2	0	3839	1	877	0	0	0
1771	0	0	0	0	1	9	0	13	0	0	0	0	0	1	1124	0	0	4	0
1811	5	1	0	12	1	0	77	0	439	0	0	0	0	6	4	5005	0	0	0
1911	151	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
1923	0	0	0	0	0	0	0	643	0	0	55	23	0	0	2	0	0	28	0
9212	0	0	0	0	23	0	0	0	0	0	0	0	0	0	0	0	0	0	33

2.2 Plot Confusion Matrix

# invert z idx values
z = cm[::-1]

x = classes_name
y = x[::-1].copy() # invert idx values of x

# change each element of z to type string for annotations
z_text = [[str(y) for y in x] for x in z]

# set up figure 
fig = ff.create_annotated_heatmap(z,
                                  x=x,
                                  y=y,
                                  annotation_text=z_text,
                                  colorscale='spectral',
                                  reversescale=True)

# add title
#fig.update_layout(title_text=f"Confusion Matrix - {site}, {year}")

# adjust margins to make room for yaxis title
#fig.update_layout(margin=dict(t=200, l=200))

#fig.update_xaxes(tickfont_size=20)
#fig.update_yaxes(tickfont_size=20)
#fig.update_layout(font_size=25)

# add colorbar
#fig['data'][0]['showscale'] = True

fig.show()
fig.write_html(cm_html, full_html=False)

3. Accuracy metrics

3.1 Compute Accuracy Metrics

If you decide that you are not interested in the scores of classes that were not predicted, then you can explicitly specify the classes you are interested in (which are labels that were predicted at least once).

acc_metrics_str  = classification_report(y_true,
                                         y_pred,
                                         target_names=classes_name,
                                         labels=classes_all,
                                         digits=3)

print(acc_metrics_str)

                                    precision    recall  f1-score   support

             Grassland and meadows      0.923     0.975     0.948     17918
                      Fruits trees      0.000     0.000     0.000       688
                         Vineyards      0.000     0.000     0.000        54
                            Forest      0.689     0.900     0.780       813
                             Urban      0.682     0.884     0.770       457
                       Greenhouses      0.167     0.029     0.049        69
          Leafy or stem vegetables      0.441     0.082     0.138       785
                      Winter wheat      0.771     0.897     0.829     22230
                             Maize      0.834     0.970     0.897      6535
                    Barley six-row      0.722     0.597     0.654      4996
                              Oats      0.208     0.372     0.267       242
                     Other cereals      0.266     0.110     0.156      4640
                          Rapeseed      0.935     0.712     0.809       650
                          Potatoes      0.995     0.790     0.881      4862
                              Peas      0.958     0.976     0.967      1152
                        Sugar beet      0.841     0.902     0.870      5550
                           Alfalfa      0.000     0.000     0.000       151
Flax  hemp and other similar crops      0.106     0.037     0.055       751
      Permanent rivers and streams      0.717     0.589     0.647        56

                          accuracy                          0.813     72599
                         macro avg      0.540     0.517     0.511     72599
                      weighted avg      0.779     0.813     0.789     72599

3.2 Compute Overall Accuracy

oa = accuracy_score(y_true, y_pred)
oa = round(oa*100, 2)

print(f'Overall Accuracy : {oa}%')

Overall Accuracy : 81.31%

3.3 Plot Accuracy Metrics

acc_metrics_dict = classification_report(y_true, y_pred,target_names=classes_name, output_dict=True)

am_df = pd.DataFrame.from_dict(acc_metrics_dict).round(3)

am_df = am_df.iloc[:,:-3]

#am_df = pd.concat([am_df, nb_df])

#am_df = am_df.sort_values(by='pix_count', ascending=False, axis=1)

class_name = am_df.columns.to_list()
precision = am_df.loc['precision'].to_list()
recall    = am_df.loc['recall'].to_list()
f1_score  = am_df.loc['f1-score'].to_list()

fig = go.Figure(data=[
    go.Bar(name='Precision', x=class_name, y=precision, text=precision, textposition='auto'),
    go.Bar(name='Recall', x=class_name, y=recall, text=recall, textposition='auto'),
    go.Bar(name='F1-score', x=class_name, y=f1_score, text=f1_score, textposition='auto')
])
# Change the bar mode
fig.update_layout(title_text=f'Accuracy Metrics - {site}, {year}',
                  barmode='group')

#fig.update_xaxes(tickfont_size=30)
fig.update_yaxes(tickfont_size=10, range=[0,1])
#fig.update_layout(xaxis_title=None, font_size=10)
#fig.update_layout(legend=dict(font=dict(size=25)))

fig.show()
fig.write_html(am_html, full_html=False)